igl, MeshCAD, Polyfem 3D example

import numpy as np
import math
import igl
import meshplot as mp
import wildmeshing as wm
import polyfempy as pf

Read the mesh

V, F = igl.read_triangle_mesh("mesh.obj")
mp.plot(V, F)

Get mesh bounding box and center

minn = np.min(V, axis=0)
maxx = np.max(V, axis=0)

center = (minn+maxx)/2

Function to select sidesets:

  • Sideset 2: z-component close to the bottom boundary minn[2] and larger than y of center
  • Sideset 3: z-component close to the bottom boundary minn[2] and smaller than y of center
def sideset(p):
    if abs(p[2] - minn[2]) < 0.5:
        if p[1] > center[1]:
            return 2
        else:
            return 3
    return 1

Wildmeshing tetrahedralization

wm.tetrahedralize("mesh.obj", "out.mesh", mute_log=True)

Loading the mesh

solver = pf.Solver()
solver.load_mesh_from_path("out.mesh")
[2019-07-22 19:36:29.077] [polyfem] [info] Loading mesh...
[2019-07-22 19:36:29.078] [geogram] [info] Loading file out.mesh...
[2019-07-22 19:36:29.140] [geogram] [info] (FP64) nb_v:4019 nb_e:0 nb_f:4872 nb_b:0 tri:1 dim:3
[2019-07-22 19:36:29.140] [geogram] [info]  nb_tets:15744
[2019-07-22 19:36:29.140] [geogram] [info] Attributes on vertices: point[3]
[2019-07-22 19:36:29.326] [polyfem] [info] mesh bb min [-92.075, -98.5272, -42.037], max [92.075, 28.6776, 42.037]
[2019-07-22 19:36:29.327] [polyfem] [info]  took 0.249492s

Using the previous function as sidesets and visualization for confirmation

solver.set_boundary_side_set_from_bary(sideset)

ps, ts, s = solver.get_boundary_sidesets()

col = np.zeros_like(s)
col[s==2] = 2
col[s==3] = 3

mp.plot(ps, ts, col, shading={"wireframe": False})

Setup the problem

settings = pf.Settings()
problem = pf.Problem()

settings.set_pde(pf.PDEs.LinearElasticity)

settings.set_material_params("E", 200)
settings.set_material_params("nu", 0.35)


problem.set_displacement(2, [0, 0, 0])
problem.set_force(3, [0, 0.5, 0])

settings.set_problem(problem)

Solve!

solver.settings(settings)
solver.solve()
[2019-07-22 19:36:29.474] [polyfem] [info] simplex_count:   15744
[2019-07-22 19:36:29.474] [polyfem] [info] regular_count:   0
[2019-07-22 19:36:29.474] [polyfem] [info] regular_boundary_count:  0
[2019-07-22 19:36:29.474] [polyfem] [info] simple_singular_count:   0
[2019-07-22 19:36:29.474] [polyfem] [info] multi_singular_count:    0
[2019-07-22 19:36:29.474] [polyfem] [info] boundary_count:  0
[2019-07-22 19:36:29.474] [polyfem] [info] multi_singular_boundary_count:   0
[2019-07-22 19:36:29.474] [polyfem] [info] non_regular_count:   0
[2019-07-22 19:36:29.474] [polyfem] [info] non_regular_boundary_count:  0
[2019-07-22 19:36:29.474] [polyfem] [info] undefined_count:     0
[2019-07-22 19:36:29.474] [polyfem] [info] total count:  15744
[2019-07-22 19:36:29.475] [polyfem] [info] Building isoparametric basis...
[2019-07-22 19:36:29.528] [polyfem] [info] Computing polygonal basis...
[2019-07-22 19:36:29.528] [polyfem] [info]  took 0.000912479s
[2019-07-22 19:36:29.531] [polyfem] [info] hmin: 0.741093
[2019-07-22 19:36:29.531] [polyfem] [info] hmax: 18.3008
[2019-07-22 19:36:29.531] [polyfem] [info] havg: 5.9001
[2019-07-22 19:36:29.531] [polyfem] [info]  took 0.0526661s
[2019-07-22 19:36:29.531] [polyfem] [info] flipped elements 0
[2019-07-22 19:36:29.531] [polyfem] [info] h: 18.3008
[2019-07-22 19:36:29.531] [polyfem] [info] n bases: 4019
[2019-07-22 19:36:29.531] [polyfem] [info] n pressure bases: 0
[2019-07-22 19:36:29.531] [polyfem] [info] Assigning rhs...
[2019-07-22 19:36:29.541] [polyfem] [info]  took 0.0102129s
[2019-07-22 19:36:29.541] [polyfem] [info] Assembling stiffness mat...
[2019-07-22 19:36:29.674] [polyfem] [info]  took 0.13286s
[2019-07-22 19:36:29.674] [polyfem] [info] sparsity: 435691/145371249
[2019-07-22 19:36:29.674] [polyfem] [info] Solving LinearElasticity with
[2019-07-22 19:36:29.676] [polyfem] [info] Hypre...

Visualization of displacement and stresses

p, t, d = solver.get_sampled_solution()
m, _ = solver.get_sampled_mises_avg()
mp.plot(p+d, t, m, shading={"wireframe": False})